A mixed formulation for a modification to Darcy equation based 
on Picard linearization and numerical solutions to large-scale 

realistic problems 

K. B. Nakshatrala and D. Z. Turner 

Abstract. In this paper we consider a modification to Darcy equation by taking into account the 
dependence of viscosity on the pressure. We present a stabilized mixed formulation for the result- 
ing governing equations. Equal-order interpolation for the velocity and pressure is considered, and 
shown to be stable (which is not the case under the classical mixed formulation). The proposed 
mixed formulation is tested using a wide variety of numerical examples. The proposed formulation 
is also implemented in a parallel setting, and the performance of the formulation for large-scale 
problems is illustrated using a representative problem. Two practical and technologically impor- 
tant problems, one each on enhanced oil recovery and geological carbon-dioxide sequestration, are 
solved using the proposed formulation. The numerical examples show that the predictions based on 
Darcy model are qualitatively and quantitatively different from that of the predictions based on the 
modified Darcy model, which takes into account the dependence of the viscosity on the pressure. In 
particular, the numerical example on the geological carbon-dioxide sequestration shows that Darcy 
model over-predicts the leakage into an abandoned well when compared to that of the modified 
Darcy model. On the other hand, the modified Darcy model predicts higher pressures and higher 
pressure gradients near the injection well. These predictions have dire consequences in predicting 
damage and fracture zones, and designing the seal, whose integrity is crucial to the safety of a 
geological carbon-dioxide sequestration geosystem. 



1. INTRODUCTION 

Darcy equation has been successfully employed to model flow in porous media in wide areas 
of applications ranging from groundwater hydrology, petroleum engineering to food engineering. 
Henry Darcy has proposed a simple equation to model the flow of an incompressible fluid in rigid 
porous media, which is popularly referred to as Darcy equation and in some instances as "Darcy 
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law." Darcy has developed this equation empirically based on his experiments on the flow of water 
through sand beds (Also see the English translation of Darcy's work by Patricia Bobeck j2].) 

It is important to note that Darcy equation is just an approximation of the balance of linear 
momentum for the fluid flowing through a rigid porous solid (see the discussion in [3, Section 2.5]). 
The equation is valid under a number of assumptions, see references [3||4] and also [§] Introduction]. 
Darcy equation merely predicts the flux, and this prediction is not accurate at high pressures. 
Moreover, in references (5|[6] it has been advocated that Darcy equation is quite good for low flow 
rates for a fluid like water but is not the case for high flow rates and for dense fluids like mineral 
oil. 

Several extensions to Darcy equation have been developed by various researchers. Two of the 
early popular ones are by Forchheimer [T] and Brinkman [s]. Bowen has outlined various models 
for flow through porous media allowing multiple fluid components and deformation of the porous 
solid |3j. Recent extensions of Darcy equation can be found in references [9-11 . A hierarchy 



of models for flow through porous media has been presented in [4] of which the simplest model is 
Darcy equation. One important point to note, which is central to this paper, is that Darcy equation 
assumes that the coefficient of viscosity (and hence the drag coefficient) to be independent of the 
pressure in the fluid. 

1.1. Pressure dependent viscosity. Certain fluid models do involve flow equations with 
non-constant viscosity (e.g., a function of pressure or velocity), which are more realistic in modeling 
enhanced oil recovery and geological carbon-dioxide sequestration. In both enhanced oil recovery 
and geological carbon-dioxide sequestration, pressure can vary from 0.1 MPa to 100 MPa. There 
is irrefutable experimental evidence that viscosity of mineral oils is not constant and changes 



drastically with respect to pressure (for example; Bridgman 12 , Andrade |13|). In fact, Barus [14| 



suggested the following exponential relationship between viscosity and pressure: 

/x(p) = /ioexp[/3p] (1) 

where /3 has units Pa^^. 

In geological carbon-dioxide sequestration, supercritical carbon-dioxide is pumped into deep 



saline aquifers or abandoned oil wells. In references 15 16 , transport properties (including vis- 
cosity) of supercritical carbon-dioxide are presented for various pressures ranging from 0.1 MPa to 
8 GPa, and for various temperatures. These experimental studies clearly indicate that the viscos- 
ity varies exponentially with respect to pressure even for supercritical carbon-dioxide for various 
temperatures and for a wide range of pressures (varying from 0.1 MPa to 8 GPa). 

It should be emphasized that petroleum reservoir simulations and geological carbon-dioxide 
sequestration involve heat transfer, mass transfer and phase changes in addition to flow aspects. 
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Specifically, one has take into account the attendant geochemical reactions (e.g., see Lichtner [17[ ). 
Several mathematical models, numerical formulations, and numerical simulations have been re- 
ported in the literature that have taken into account more than the flow aspects. For example, see 



references [18-20 . The mathematical model in this paper does not address such issues but just 
focuses on the need to take into account the dependence of viscosity on the pressure in enhanced 
oil recovery and geological carbon-dioxide sequestration simulations. Our study is primarily moti- 
vated by the fact that the flow characteristics of a fluid whose viscosity depends on pressure can be 
significantly different from that of the flow characteristics of a fluid with constant viscosity. Some 
representative prior works on modeling of fluids by taking into account the dependence of viscosity 
on pressure are 21-24 . In this paper we consider a modification to Darcy equation that takes 
into account the dependence of viscosity on pressure. The resulting equations will he in mixed form 
and nonlinear. The unknowns are the velocity and pressure. We present a new stabilized mixed 
formulation based on the variational multiscale formalism and fixed-point linearization. 

1.2. Stabilized mixed finite element formulations. It is well-known that care should 
be taken to avoid numerical instabilities when dealing with mixed formulations. As discussed by 
Franca and Hughes |25|, a stable mixed formulation either meets or circumvents the Ladyzhenskaya- 



Babuska-Brezzi (LBB) inp-sup stability condition [26,27 . Stabilized methods typically fall in the 
later category, and the stability is achieved through addition of stabilization terms. But most 
of the approaches that aim to satisfy LBB condition achieve the stability by restricting the in- 
terpolation functions for the independent variables. For example, PlPO approach in which the 
velocity unknowns are placed at nodes, and pressure unknowns are for each element /cell. Two 
other notable works to satisfy the LBB condition are to use either the Raviart-Thomas spaces [28| 
or Brezzi-Douglas-Marini spaces [29] . In our paper, we do not take such an approach of placing 
restrictions on the interpolation functions for the independent variables. Instead, we circumvent 
the LBB condition by adding stabilization terms in a consistent manner, and allow any combination 
of interpolation functions for the independent variables (in this case, the velocity vector and the 
pressure) . 

It should be emphasized that the LBB stability condition is applicable even to mixed formu- 
lations based on the finite volume method. As mentioned earlier, a mixed method has to either 
satisfy the LBB condition in the case of saddle-type formulation, or circumvent the LBB condition 
by adding stabilization terms to avoid saddle-type problem. In the context of a mixed formula- 
tion for Darcy-type equation, the LBB stability condition places a restriction on the choice of the 
function spaces for the vector-field variable (i.e., the velocity) and the scalar-field variable (i.e., the 
pressure). It is interesting to note that staggered grids, which is quite popular in the finite volume 
literature, is a way to satisfy the LBB condition. 



Some popular approaches for developing stabilized finite element formulations are least-squares 



|30|[3l] , Galer kin/least-squares (GLS) [32t[33] , finite increment calculus (FIC) [34t[35] , streamline 
upwind Petrov-Galerkin |36|. Another popular approach for developing stabilized formulations is 
the variational multiscale framework |37|j, which has been employed in this paper. The variational 
multiscale framework has been successfully employed in many studies to develop stabilized formu- 
lations for a wide variety of problems: Darcy equation 38-40 , Stokes' equation (41 42 , linearized 
elasticity (43,44 , incompressible Navier-Stokes 45,46 , Fokker-Planck |47| . 

A huge volume of literature is available on mixed methods and stabilized formulations, and a 
thorough discussion on these topics is beyond the scope of this paper. Some representative papers 



on stabilized formulations are 125 , 32 , 36 , 38 , 40 , 48 -551. Some texts concerning stabilized methods 



are (56-59 . Some representative works on mixed methods in the context of flow through porous 
media are (38H40l[60H67] . It should be noted that none of the aforementioned numerical works 
considered the fluid model considered in this paper. 

1.3. Main contributions of this paper. To the best of our knowledge no other prior stud- 
ies have systematically studied (in a numerical setting) the modified Darcy equation except for 
Reference [s] . But this paper presents an alternate and simpler formulation than the one presented 
in Reference (s]. The simplicity arises due to the fact that the proposed formulation is obtained 
by first linearizing the governing equations and then applying the variational multiscale formalism 
on the resulting equations. This approach facilitates to extend the proposed formulation to other 
complicated models (e.g., modified Brinkman by taking into account the pressure dependent vis- 
cosity, Brinkman-Forchheimer) in a straightforward manner, which may not be the case with the 
stabilized mixed formulation proposed in Reference [S]). The main contributions of this paper can 
be summarized as follows: 

(a) We considered a realistic modification to Darcy equation by taking into account the dependence 
of viscosity on the pressure. This implies that the drag coefficient will depend on the pressure. 
We presented a new stabilized mixed formulation based on the variational multiscale formalism 
for the resulting nonlinear equations. 

(b) We have shown numerically that equal- order interpolation for the velocity and pressure (which is 
computationally the most convenient) is stable under the proposed formulation. A wide variety 
of test problems are performed to illustrate the performance of the proposed formulation. 

(c) We have implemented the proposed in a parallel setting. 

(d) We have solved two practical and technologically important large-scale problems with relevance 
to enhanced oil recovery and geological carbon-dioxide sequestration. The numerical results 
have clearly indicated the importance of considering the role of dependence of viscosity on the 
pressure in these two application areas. 



1.4. An outline of the paper. The remainder of the paper is organized as follows. In 
Section [2| we present the modified Darcy equation, and present the proposed stabilized mixed 
finite element formulation, which is based on the variational multiscale formalism. We also present 
a numerical solution procedure for solving the resulting nonlinear equations. In Section [Sj we 
illustrate the performance of the proposed formulation on a wide variety of benchmark problems, 
which are commonly used in the literature for testing mixed formulations. In the same section 
we also present the numerical results for two large-scale practical problems, which are solved by 
implementing the proposed formulation in a parallel setting. Conclusions are drawn in Section [4| 

2. GOVERNING EQUATIONS AND STABILIZED FORMULATION 

Let O C M"°' be an open and bounded set, where "nd" denotes the number of spatial dimensions. 
Let d^l := Cl — Qhe the boundary (where 0, is the set closure of fi), which is assumed to be piecewise 
smooth. A spatial point in 0, is denoted by x. The spatial gradient and divergence operators are, 
respectively, denoted as "grad[-]" and "div[-]". Let v : Q, ^ M""^ denote the velocity field, and 
p : — )■ M denote the pressure field. The boundary is divided into two parts, denoted by 
and r^, such that F'^ n F^' = and F'' U F^ = dQ. is the part of the boundary on which 
normal component of the velocity is prescribed, and F^ is part of the boundary on which pressure 
is prescribed. 

The modified Darcy equations can be written as 

a{p)v + grad[p] = p{x)b{x) in Q. (2a) 

div[^;] = in J] (2b) 

v{x) • n{x) = Vn{x) on F" (2c) 

p{x) = po{x) on TP (2d) 

where a{p) is the drag coefficient (which has dimensions of [ML^'^T^^]), po{x) is the prescribed 
pressure, Vnix) is the prescribed normal component of the velocity, p{x) denotes the density of the 
fluid, b{x) is the specific body force, and n{x) is the unit outward normal vector to dO,. The drag 
function is the ratio between viscosity of the fluid and the permeability. Herein we consider the 
following two forms for drag function 

a{p) = ao(l + Pp) (3a) 
a{p) = aoexp[f3p] (3b) 

The first equation can be considered as a two-term Taylor's series approximation of the second 
equation (which is based on the Barus' formula). It is, in general, not possible to obtain analytical 
solutions for the system of equations Q, especially for complex geometries. Hence, one may have 
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to resort to numerical solutions. The main aim of this paper is to present a stabilized mixed 



formulation to solve the boundary value problem given by equations (2a)-(2d). 

Developing numerical formulations for the above equation falls in the realm of mixed meth- 
ods [27| . As mentioned earlier, it is generally agreed upon that care should be taken in developing 
numerical formulations to avoid numerical instabilities. For example, under the classical Galerkin 
formulation (which is sometimes referred to as the classical mixed formulation) equal-order inter- 
polation for velocity and pressure is not stable. A mathematical theory that addresses the stability 
issues with mixed methods is the Ladyzhenskaya-Babuska-Brezzi (LBB) inf-sup condition. Special- 
ized elements and appropriate function spaces have been developed for mixed problems (like Darcy 
equation, Stokes equation, incompressible Navier-Stokes) to either satisfy or circumvent the LBB 
stability condition. 

In the next subsection we present a new stabilized mixed formulation that is inherently more 
stable than the classical mixed formulation and can accommodate a larger variety of combinations of 
interpolation functions for velocity and pressure including the equal-order interpolation. Although 
the proposed formulation can accommodate many combinations of interpolations for the velocity 
and pressure, in Section [3] we illustrate the performance of the proposed formulation using equal- 
order interpolation as other combinations are not computationally attractive. 

2.1. A new stabilized mixed formulation for modified Darcy equation. Let w{x) and 
q{x) be test functions corresponding to the velocity and pressure, respectively. Let us define the 
following function spaces, which will be used in the remainder of this paper: 

V := \^v{x) G {L2{n))'"^ I d[v[v] G L2{n),v{x) ■ n{x) = Vn{x) on (4a) 
W := \^w{x) e {L2{n))""^ I divH G L2in),w{x) • n{x) = on (4b) 
Q = H\n) (4c) 



In using equation (4c) one assumes that 7^ 0. If F-^ = then, for well-posedness, Q should be 
taken as follows: 



Q ■= |p(a;) G H\n) \ j p{x) df] = o| 



In addition, for the case of F^ = 0, we need to satisfy the following compatibility condition on the 
prescribed data: 



/ Vn{x)dT = {) (5) 

Jan 



Ian 
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For convenience, let us denote the standard L2 inner product defined over spatial domain K as 
(•; ■)k- That is, 

{a;b)K= [ a-bdK ya,b (6) 
Jk 

For simplicity, the subscript K will be dropped if K = Q. 



In References 38,40| a stabilized formulation based on the variational multiscale (VMS) for- 
malism has been proposed and analyzed for the case of Darcy equation (which assumes constant 
drag coefficient). In this paper we extend the variational multiscale formulation to the modified 
Darcy equation (which gives rise to nonlinear equations) . We shall linearize the governing equations 



(2a)-(2d) using a fixed-point procedure as follows: 

+ grad[p(^)] = p{x)b{x) in n (7a) 

dw[v^'^] =0 in 17 (7b) 

v'^'\x) ■ n{x) = Vn{x) inP'' (7c) 

p('\x)=po{x) on TP (7d) 



Based on the derivation presented in Reference |40| , a stabilized mixed formulation for the above 
equations based on the VMS formalism can be written as follows: Find v^^\x) G V and p^^\x) G Q 
such that we have 

gstMw,q;v^'\p^'^;p^'-'^) = Cst.b{w,q;p^'~'^) yw{x) G W,q{x) G Q (8) 
where the functionals ^stab and £stab are, respectively, defined as follows: 
Gst!ih{w,q;v,p;p) := {w;a{p)v) - {div[w];p) - {q;dw[v]) 

- ^ {oi{p)w + grad[q']; a"-'^(p) {a{p)v + grad[p])) (9a) 
£-stah{w,q;p) := {w; pb) - {w ■ n;po)rp - ^ {a{p)w + gTad[q];a''^{p)pb) (9b) 



The terms in the second line of equation ( |9a[ ) and the last term in (9b) are referred to as stabilization 
terms. The factor 1/2 (in front of these terms) is the stabilization parameter. The last term in 



equation (9b) is the stabilization term due to the body force. The stability of the above formulation 



can be inferred using the mathematical proof outlined in Reference 38 



The proposed numerical algorithm for solving the original system of equations ([2]) is given in 

Algorithm [T| and will be based on the above stabilized formulation. Since the formulation is based 

on the fixed-point linearization, the rate of convergence of the algorithm (if it converges) will be 

linear with respect to iteration number. Moreover, when a is a constant (that is, for classical 

Darcy equation), the proposed formulation will converge in one iteration. In a finite element 
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implementation, the norm in the stopping criterion (see hne 10 in Algorithm [T]) can be taken as 
the Euclidean norm of the vector containing nodal pressures. 

Algorithm 1 A new stabilized formulation for modified Darcy equation 
1: Input: eTOL, MAXITERS 
2: Output: Vh{x) and ph{x) 
3: Guess p^l^\x) € Qh 
4: Initialize: i 1 
5: while true do 
6: if i > MAXITERS then 

7: print Iterative scheme did not converge. RETURN 
8: end if 

9: Find v^^^^ (x) G Vh and pj^*"* (a;) G Qh such that 

Qsta.h{w h,qh;v''f^\ph^;Ph'^^) = Csti,h{wh,qh;Ph~^^) ^Wh{x) g Wh, qh{x) g Qh 

{Check for convergence based on the given tolerance} 
10: if \\pf -Ph'^^'W < exoL then 
11: print Iterative scheme converged. 
12: Vh{x) ^ v^f^\x) and ph{x) ^ Ph\x) 
13: RETURN 
14: end if 
15: i ^ i + 1 
16: end while 



3. REPRESENTATIVE NUMERICAL RESULTS 

In this section, we present a wide variety of representative numerical results to illustrate the 
performance of the proposed stabilized formulation. The first few test problems are used to illus- 
trate the accuracy of the proposed formulation on various computational grids as for these test 
problems analytical solutions can be easily obtained. One test problem is used to illustrate the 
robustness of the formulation in the case of heterogeneous medium. To the end of the section, 
the performance of the proposed formulation in a parallel setting is illustrated using a large-scale 
simulation related to carbon-dioxide sequestration. It should be noted that in all our numerical 
simulations we have employed low-order Lagrange finite element approximations (i.e., p= 1), and 
equal-order interpolation for the velocity and pressure. 

We shall first non-dimensionalize the governing equations by employing the same non-dimensionalization 
procedure as given in Reference [s]. We shall present the details for completeness. We define the 
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following non-dimensional quantities (which have a superposed bar): 

X _ V _ p _ a _ ao _ P T b Q OTD - PO/inA 

= y, I' = TT, P = " = ' "0 = ' P = ' ^ = P = '"n = jy, Po = s (10) 

L V F Uj-cl Oref Prcf B VP 

where L, V , P, Qref, Pref and B respectively denote reference length, velocity, pressure, drag 
coefficient, density and specific body force. The gradient and divergence operators with respect to 
X are denoted as "grad" and "div", respectively. The scaled domain ilscaicd is defined as follows: a 
point in space with position vector x G f^scaied corresponds to the same point with position vector 
given by a; = xL G Q . Similarly, one can define the scaled boundaries: ^fiscaied, r^^^^g^, and T^^^jg^. 
The above non-dimensionalization gives rise to two dimensionless parameters 

A:=^^ and C := ^ (11) 

A corresponding non-dimensionalized form of the drag functions given in equation ([s]) can be written 

as 



a{p) = ao (1 + /3p) (12a) 
q(p) = Qoexp[/3p] (12b) 



We shall write a non-dimensional form of the modified Darcy equation as 



A a{p)v + grad[p] = C p h{x) in l^scaied (13a) 

div[^)] = in riscaicd (13b) 

v{x) ■ n{x) = Vn{x) onT^g^igd (13c) 

P{x)=Pq{x) onr°^igd (13d) 

where fi{x) is the unit outward normal to the boundary f^Ogcaled- 

3.1. One-dimensional problem. We shall test the proposed stabilized mixed formulation 
using a simple one-dimensional problem. Consider the computational domain to be of unit length, 
and pressures of pi and p2 are respectively prescribed at the left and right ends of the unit domain 
(see Figure [1]). We neglect the body force. The governing equations for this test problem can be 
written as 

Aa{p)v{x) ^ = in (0, 1) (14a) 
ax 

du 

— = in (0, 1) (14b) 
ax 

p{x = ff)=pi, p{x = l)=p2 (14c) 
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For this test problem, the analytical solutions with the drag functions defined in equation ([3| can 
be written as 



for the case a{p) = ao{l + /3p) < 



/9 

v(x) = -j^r^ In 



l+/3p2 



:i5a) 



p(a^) = -f" In { (1 — x) exp [— + X exp [— /3p2] } 
for the case q:(p) = ao exp[/3p] < _ _ (15b) 

In Figures [2] and [3] we compare the numerical solutions against the analytical solutions for various 
drag functions. In all the cases considered, the proposed numerical formulation performed well. It 
is important to note that steep gradients in the pressure occur near the boundary for high values of 
/3 (which is an indicator of the strength of nonlinearity) . Figure |4] shows the number of iterations 
taken by the numerical algorithm with respect to /?. As expected, the number of iterations increases 
with respect to (3. Figure |4] shows the variation of the residual with respect to the iteration number 
for various value of (3. The residual decreases monotonically with the iteration number. 

3.2. Two-dimensional constant flow problem. The computational domain is a bi-unit 
square domain, = [0, 1] x [0, 1]. Barus formula with ao = 1 is employed with two different values 
of /3 = 0.1 and f3 = 0.4. The body force is neglected. The left and right sides of the domain are 
prescribed with Vx = 1- The top and bottom sides of the domain are prescribed with vg = 0. The 
analytical solution is given by 

Vx{x,y) = l, Vy{x,y) = (16a) 
p{x, y) = po - i In [1 - ao^(l - x) exp[^po]] (16b) 

We prescribed po = 1 at the top-right corner of the computational domain. The finite element 
meshes used in the numerical simulations are shown in Figure |6j The pressure contours using the 
proposed stabilized formulation are shown in Figure [TJ 

3.3. The quarter five-spot problem. A standard test problem is the quarter five-spot prob- 
lem. The computational domain is a square in the horizontal plane with injection and production 
wells at opposite corners along one of the diagonals. The source and sink strengths at injection and 
production wells are, respectively, taken as -1-1/4 and —1/4. There is no volumetric source/sink (i.e., 
b{x) = 0). The five-spot problem is pictorially described in Figure [Sj Figure [9] show the pressure 
contours using the proposed stabilized formulation. In these simulations we have employed Barus' 
formula (12b) with cxq = 1, ^ = 0.3, A = 1, and B = 0. The variation of the norm ~ Ph ^^11 



with respect to iteration number is shown in Figure [10} As one can see from these figures, there are 



no spurious oscillations in the pressure, and the proposed stabilized formulation performed well. 
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Figure |4] shows that the number of iterations increase with an increase in /3, which is a measure 
of the strength of the nonhnearity. Figure |4] shows that the norm of the residual vector decreases 
monotonically with respect to the iteration number under the proposed stabiUzed formulation for 
the five-spot problem. 

3.4. The checkerboard problem. This problem tests the formulation for the case in which 
there are abrupt changes in the drag coefficient. The geometry and boundary conditions are same 
as the quarter five-spot test problem. But, for this test problem, the computational domain is 



divided into four regions as shown in Figure 13 In regions I and IV, we have taken qq = 1; and 
in regions II and III, we have taken ao = 0.001. In the numerical simulations we have employed 



Barus formula with f3 = 0.3, A = 1, B = 0, and have taken etol = 10 ^. In Figure 14 we have 
shown the contours of pressure for both four-node quadrilateral and three-node triangular meshes. 



-(i) -(i~l) 

In Figure 15, we have plotted the variation of the norm ~ Ph II '^ifh respect to iteration 



number. As one can see from these figures, there are no spurious oscillations in the pressure field, 
and the proposed stabilized formulation performed well. 

3.5. Numerical /i-convergence studies. Consider the computational domain to be = 
(0, 1) X (0, 1). We have employed Barus formula with ao = 1 ^-nd f3 = 0.5. The exact solution for 
the velocity and pressure is given by 

Vx{x, y) = + sin(7rx) cos(7ry), Vy{x, y) = - cos{-kx) sin(7ry), p{x, y) = x^y"^ (17) 

The density is taken to be unity, and the specific body force is given by 

hx{x^ y) = + exp(x^y^/2) sin(7rx) cos(7ry) + 2xy^, by{x, y) = — exp(x^y^/2) cos(7rx) sin(7ry) -|- 2x^y 

(18) 

The normal component of the velocity is prescribed to be zero on the boundary, and etol = 10^^. 



The rate of convergence with respect to mesh refinement is shown in Figure 16 The pressure 



contours are shown in Figure 17 From these figures it is evident that the proposed algorithm 
performed well. 

3.6. Three-dimensional constant flow. This problem is a three-dimensional version of the 
patch test, which is employed in Reference [40] to assess the stability of a stabilized formulation 
for Darcy equation. Herein, we shall use the patch test to assess the stability of proposed stabilized 
formulation for modified Darcy equation. The formulation is considered to pass the patch test if 
the flow field matched with the analytical solution up to machine precision. 

The computational domain is a cube given by = (0,5) x (0,5) x (0,5). We have employed 

Barus formula with = 1 and /3 = 0.1. We have prescribed the normal component of the velocity 

on x = and x = 5 faces to be unity. On the other four outer faces we have prescribed f„ = 0. 
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The body force is assumed to be zero, and for uniqueness the pressure at the origin is taken to be 
zero (that is, po(0, 0, 0) = 0). The analytical solution for the pressure is given by 

p{x,y,z) = -^ln[l + aof3x] (19) 



We have taken exoL = 10^^ in this numerical simulation. In Figure 18 we have shown the pres- 
sure contours obtained using the proposed formulation, and the numerical matched well with the 
analytical solution. 

3.7. Regions with different permeability. This problem considers fully saturated, single 
phase, single component flow in the region of a production well close to a boundary between regions 
of permeability that differ by orders of magnitude. Figure ?? shows the problem domain, the 
production well near the straight interface, and the computational mesh employed in the numerical 
simulation. The mesh consists of 12, 924 eight-node brick linear elements with a total of 63,224 
unknowns. The smallest elements (which are near the well) are 0.33 units in diameter while the 
largest elements (which are near the perimeter of the circle) are 5.0 units in diameter. We have 
taken eTOL = 10^^^. The properties used in this analysis are listed in Table [ij The boundary 
conditions for this problem consists of no-penetration boundaries at the top and bottom of the 
domain (except for the opening at the production well). The pressure at the well opening and 
around the sides is prescribed in a weak fashion with values given in Table [T] The drag function 
a{p) is calculated using equation (12a) with /3 varying between 0.0 and 1.0. 

The exact solution for the velocity streamlines for this problem for constant drag function can be 



obtained by the method of images and is presented in Reference 68 . In region A, the streamlines 
should curve away from the production well in a radial fashion. In region B, the streamlines should 



remain straight. Figure 21 shows similar results for the streamlines as produced by the method 
presented in this work. Note that the streamlines in Region B are straight as predicted by the 



exact solution. Figure 22 shows a log-log plot of the production rate at the well as a function of 
f3. The production rate was computed as the total flux of fluid J v{x) ■ n{x) dF across the well 
opening. Notice that as f3 increases, the production rate decreases in a linear fashion. Figure 23 
shows the amount of the total production at the well that emanates from region A or B. This result 
is of particular interest in determining the effect of {3 on the amount which the well will draw from 
either region. The results show that although as f3 increases the total production decreases, the 
proportion of the total production that comes from region A or B remains constant. 

3.8. CO2 leakage through an abandoned well (large-scale problem). The last numer- 
ical example makes an important contribution to the study of geological carbon-dioxide (CO2) 
sequestration into underground aquifers. When CO2 is pumped into an underground aquifer, it 

can leak through fissures in the surrounding aquitard layers or though man-made penetrations, 
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Table 1. Parameters used in regions with different permeability problem 



Parameter Value 

ao region A 0.001 

OfQ region B 1.0 

p 1.0 

6 (0,0,0) 

Po far- field boundary 1.0 

Po well opening 3.333 xlO^^ 

Radius of far-field boundary 100.0 

Radius of well 1.0 



such as abandoned wells. An important question regarding the suitability of certain locations 
for sequestration involves predicting the leakage rate of CO2 into other aquifers. This numerical 
example models the leakage of CO2 through an abandoned well as it is injected into an aquifer. 



The geometry is shown in Figure 24, The diameter of the wells are 0.3 m and they are located 
100 m apart from each other. The far-field boundary is located at a radius of 500 m. The com- 



putational mesh for this problem, which shown in Figure 25, consists of 271,050 eight-node brick 
linear elements and 284,019 nodes for a total of 1.14 million unknowns. Near the injection well 
and the abandoned well the element diameter is 0.05 m. Near the far-field boundary the element 
diameter is 2.5 m. To simplify the problem we assume fully saturated, single component, single 
phase incompressible fiow and that there is no body force, b = (0,0,0). The boundary conditions 
consist of no penetration boundaries on the top and bottom of the aquifers and weakly prescribed 
pressure boundary conditions along the sides. For aquifer B, a constant reference pressure, pref, 
of 2.9315 X 10^ Pa is prescribed along the outer boundary. For aquifer A, a constant reference 
pressure of 3.0599 x 10'' Pa is prescribed. At the injection well, a constant infiow velocity of 0.262 
m/s is prescribed. The parameters used for this problem are listed in Table [2j a{p) is determined 
according to the exponential function given in equation (12b), with /3 varying from 0.0 to 1.0. We 
have again taken etol = 10~^^. 

This problem was solved in parallel setting with 32 processors, on the Tri-Lab Linux Capacity 



Cluster, using Aria computer code 69 . Load balancing across the processors was achieved using 



the Zoltan [70| package. The total CPU time for a single run with a given (3 was 35 minutes with 
an average memory use on each processor of 400 MB. A complete investigation into the parallel 
scalability of this algorithm is intended for future work. 



Figure 26 shows contours of the pressure for the minimum and maximum values of j3. Notice 



that as /3 increases, the pressure in the region surrounding the injection well increases. Figure 27 
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shows the magnitude of the velocity for the maximum and minimum values of /?. As /3 increases, 
the velocity magnitude in the region of the abandoned well decreases suggesting that not including 



the pressure dependent viscosity under predicts the leakage rate. Figure 28 shows the ratio of 
injection rate to leakage rate as /3 increases. The injection rate was computed as the total flux of 
fluid J v{x) ■n(x) dT across the injection well opening. Likewise, the leakage rate was computed as 
the total flux of fluid across a cross-section of the abandoned well at the bottom of the aquifer B. 
In References |71[|72| , the authors present transient results for a similar simulation for multiphase 
flow of CO2 and brine, with different saturations, and constant viscosity. They observe peak ratios 
of leakage rate to injection rate in the range of 0.2% to 0.42%. For the simulation presented in this 
work, the ratio is between 1% and 10%. Clearly, the assumptions of fully saturated, single phase 
flow do not accurately capture the leakage rate, but the results are meaningful in that they show, 
in general, an over- prediction of the leakage rate if a pressure dependent viscosity is not accounted 
for. As /3 increases, the amount CO2 that leaks into aquifer B is decreased significantly. Such 
results suggest that a pressure dependent viscosity model has a substantial effect on the predicted 
leakage rate. 

Table 2. CO2 leakage through an abandoned well: problem parameters 



Parameter Value 

ao aquifer A 1.0 

ao aquifer B 1.0 

ao inside the wells 100 

Po aquifer A 1.0 

po aquifer B 0.95 

13 1.0 X lO"'^ 

V ■ n at inflow 0.262 



4. CONCLUDING REMARKS 

We have considered a modification to Darcy equation by taking into account the dependence 
of viscosity on the pressure, which has been observed in many experiments on organic liquids. We 
have developed a new mixed stabilized formulation for the modified Darcy equation. We have also 
presented a numerical solution procedure to solve the resulting nonlinear equations. It has been 
shown numerically that equal-order interpolation for the velocity and pressure is stable under the 
proposed stabilized mixed formulation, which is not the case with the classical mixed formulation. 
Using representative problems, it has been shown that 
(a) the proposed stabilized mixed formulation performs well, and 
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(b) the results predicted by the standard Darcy model are qualitatively and quantitatively different 
from that of the predictions based on the modified Darcy model. 

It has been observed that the dependence of viscosity on the pressure drastically alters the pressure 
profile in the domain, and creates steep gradients near the boundary of the domain. Using a repre- 
sentative problem with relevance to geological carbon-dioxide sequestration, it has been shown that 
the standard Darcy model over-predicts the velocity of the fluid in the abandon well in comparison 
with the modified Darcy model. This prediction will serious consequences in designing the ceiling 
of the cap rock, which is one of the main mechanisms for the safety of a geological carbon-dioxide 
geosystem. Another important point to be noted is that the modified Darcy model predicts higher 
pressures than that of the Darcy model, which will have implications in modeling damage and 
fracture of the porous solid. 

As a part of future work on the numerical front, one can extend the current work to solve 
problems with much higher values by employing continuation-type methods (that is, to solve the 
problem for a high fi using information from the solution at a lower fi) . Another interesting future 
numerical work can be designing preconditioners for these kinds of nonlinear problems. A possible 
future work on the modeling front is to model the damage and fracture of the porous solid along 
with the flow aspects by taking into account the dependence of the viscosity on the pressure. 

ACKNO WLED GMENTS 

The first author (K. B. Nakshatrala) was supported in part by the Department of Energy 
through the Subsurface Biogeochemical Research (SBR) Program. The authors also acknowledge 
the financial support from Sandia National Laboratories through the Laboratory Directed Research 
and Development program. The opinions expressed in this paper are those of the authors and do 
not necessarily reflect that of the sponsors. The authors thank Pat Notz and Mario Martinez for 
their insightful discussions on this topic. 

References 

[1] H. Darcy. Les Fontaines Publiques de la Ville de Dijon. Victor Dalmont, Paris, 1856. 

[2] P. Bobeck. The Public Fountains of the City of Dijon. Kendall Hunt Publishing, Dubuque, Iowa, USA, 2004. 
[3] R. Bowen. Porous Elasticity: Lectures on the Elasticity of Porous Materials as an Application of the Theory of 

Mixtures. Available online: http://rcpository.tamu.edu/handle/1969. 1/2500, Texas A&M University, 2010. 
[4] K. R. Rajagopal. On a hierarchy of approximate models for flows of incompressible fluids through porous solids. 

Mathematical Models and Methods in Applied Sciences, 17:215-252, 2007. 
[5] K. B. Nakshatrala and K. R. Rajagopal. A numerical study of fluids with pressure dependent viscosity flowing 

through a rigid porous medium. International Journal for Numerical Methods in Fluids, DOI: 10.1002/fld.2358, 

2010. 

15 



[6] D. Munaf, A. S. Wineman, K. R. Hajagopal, and D. W. Lee. A boundary value problem in groundwater motion 
analysis-comparisons based on Darcy's law and the continuum theory of mixtures. Mathematical Modeling and 

Methods in Applied Science, 3:231-248, 1993. 
[7] P. Forchheimer. Wasserbewegung durch Boden. Zeitschrift des Vereines Deutscher Ingeneieure, 45:1782-1788, 
1901. 

[8] H. C. Brinkman. On the permeability of the media consisting of closely packed porous particles. Applied Scientific 
Research, Al:81-86, 1947. 

[9] J. M. D. Medeiros, J. M. Gurgel, and F. Marcondes. Numerical analysis of natural convection in porous me- 
dia: The influence of non-Darcian terms and thermal dispersion. Journal of Porous Media, 9:235-250, 2006. 
[10] K. Al-Farhany and A. Turan. Non-Darcy effects on conjugate double-diffusive natural convection in a variable 
porous layer sandwiched by finite thickness walls. International Journal of Heat and Mass Transfer, 54:2868- 
2879, 2011. 

[11] T. R. Mahapatra, D. Pal, and S. Mondal. Influence of thermal radiation on non-Darcian natural convection in a 
square cavity filled with fluid saturated porous medium of uniform porosity. Nonlinear Analysis: Modelling and 
Control, 17:223-237, 2012. 

[12] P. W. Bridgman. The Physics of High Pressure. MacMillan Company, New York, USA, 1931. 

[13] E. N. da. C. Andrade. Viscosity of liquids. Proceedings of the Royal Society of London, Series A, 215:36-43, 
1952. 

[14] C. Barus. Isotherms, isopiestics and isometrics relative to viscosity. American Journal of Science, 45:87-96, 1893. 
[15] V. Vesovic, W. A. Wakeham, G. A. Olchowy, J. V. Sengers, J. T. R. Watson, and J. Millat. The transport 

properties of carbon dioxide. Journal of Physical and Chemical Reference Data, 19:763-808, 1990. 
[16] E. H. Abramson. Viscosity of carbon dioxide measured to a pressure of 8 GPa and temperature of 673 K. Physical 

Review E, 80(021201), 2009. 

[17] P. C. Lichtncr. Continuum formulation of multicomponent-multiphase reactive transport. Reviews in Mineralogy 

and Geochemistry, 34:1-81, 1996. 
[18] P. C. Lichtner. Continuum model for simultaneous chemical reactions and mass transport in hydrothermal 

systems. Geochimica et Gosmochimica Acta, 49:779-800, 1985. 
[19] C. Lu and P. C. Lichtner. High resolution numerical investigation on the effect of convective instability on long 

term CO2 storage in saline aquiferes. Journal of Physics: Conference Series, 78:012042, 2007. 
[20] F. Marcondes and K. Scpehrnoori. An clement-based finite-volume method approach for heterogeneous and 

anisotropic compositional reservoir simulation. Journal of Petroleum Science and Engineering, 73:99-106, 2010. 
[21] J. Malek, J. Necas, and K. R. Rajagopal. Global analysis of the flows of fluids with pressure-dependent viscosities. 

Archive for Rational Mechanics and Analysis, 165:243-269, 2002. 
[22] S. C. Subramaniam and K. R. Rajagopal. A note on the flow through porous solids at high pressures. Computers 

and Mathematics with Applications, 53:260-275, 2007. 
[23] M. Bulicek, J. Malek, and K. R. Rajagopal. Navier's slip and evolutionary Navier-Stokes-like systems with 

pressure and shear-rate depedent viscosity. Indiana University Mathematics Journal, 56:51-85, 2007. 
[24] K. Kannan and K. R. Rajagopal. Flow through porous media due to high pressure gradients. Applied Mathematics 

and Computation, 199:748-759, 2008. 
[25] L. P. Franca and T. J. R. Hughes. Two classes of mixed flnite element methods. Computer Methods in Applied 

Mechanics and Engineering, 69:89-129, 1988. 



16 



I. Babuska. Error bounds for finite element methods. Numerische Mathematik, 16:322-333, 1971. 

F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods, volume 15 of Springer series in computa- 
tional mathematics. Springer- Verlag, New York, USA, 1991. 

P. A. Raviart and J. M. Thomas. A mixed finite element method for 2nd order elliptic problems. In Mathematical 
Aspects of the Finite Element Method., pages 292-315, Springer- Verlag, New York, 1977. 

F. Brezzi, J. Douglas, and L. D. Marini. Two families of mixed elements for second order elliptic problems. 
Numerische Mathematik, 47:217-235, 1985. 

P. B. Bochev and M. D. Gunzburger. Least-Squares Finite Element Methods. Springer, New York, USA, 2009. 

B. Jiang. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics 

and Electromagnetics. Springer- Verlag, New York, USA, 1998. 

T. J. R. Hughes, L. Franca, and G. Hulbert. A new finite element formulation for computational fluid dynam- 
ics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Computer Methods in Applied 
Mechanics and Engineering, 73:173-189, 1989. 

C. Baiocchi, F. Brezzi, and L. Franca. Virtual bubbles and Galerkin-least-squares type methods (Ga.L.S.). 
Computer Methods in Applied Mechanics and Engineering, 105:125-141, 1993. 

E. Onate. A stabilized finite element method for incompressible viscous fiows using a finite increment calculus 
formulation. Computer Methods m Applied Mechanics and Engineering, 182:355-370, 2000. 
E. Onate. Possibilities of finite calculus in computational mechanics. International Journal for Numerical Methods 
in Engineering, 60:255-281, 2004. 

A. N. Brooks and T. J. R. Hughes. Streamline-upwind/Petrov-Galerkin methods for convection dominated flows 
with emphasis on the incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and 
Engineering, 32:199-259, 1982. 

T. J. R. Hughes. Multiscale phenomena: Green's functions, the Dirichlct-to-Neumann formulation, subgrid scale 
models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering, 
127:387-401, 1995. 

A. Masud and T. J. R. Hughes. A stabilized mixed flnite element method for Darcy flow. Computer Methods in 
Applied Mechanics and Engineering, 191:4341-4370, 2002. 

T. J. R. Hughes, A. Masud, and J. Wan. A stabilized mixed discontinuous Galerkin method for Darcy flow. 

Computer Methods in Applied Mechanics and Engineering, 195:3347-3381, 2006. 

K. B. Nakshatrala, D. Z. Turner, K. D. Hjelmstad, and A. Masud. A stabilized mixed finite element formulation 
for Darcy fiow based on a multiscale decomposition of the solution. Computer Methods m Applied Mechanics 
and Engineering, 195:4036-4049, 2006. 

A. Masud. A stabilized mixed finite element method for Darcy-Stokes flow. International Journal for Numerical 
Methods in Fluids, 54:665-681, 2007. 

D. Z. Turner, K. B. Nakshatrala, and K. D. Hjelmstad. On the stability of bubble functions and a stabilized 
mixed finite element formulation for the Stokes problem. International Journal for Numerical Methods in Fluids, 
60:1291-1314, 2009. 

A. Masud and K. Xia. A stabilized mixed finite element method for nearly incompressible elasticity. Journal of 
Applied Mechanics, 72:711-720, 2005. 

K. B. Nakshatrala, A. Masud, and K. D. Hjelmstad. On flnite element formulations for nearly incompressible 
linear elasticity. Computational Mechanics, 41:547-561, 2008. 



17 



[45] A. Masud and R. A. Khurram. A multiscale finite element method for the incompressible Navier-Stokes equations. 

Computer Methods in Applied Mechanics and Engineering, 195:1750-1777, 2006. 
[46] D. Z. Turner, K. B. Nakshatrala, and K. D. Hjelmstad. A variational multiscale Newton-Schur approach for 
the incompressible Navier-Stokes equations. International Journal for Numerical Methods in Fluids, 62:119-137, 
2010. 

[47] A. Masud and L. A. Bergman. Application of multiscale finite element methods to the solution of the Fokker- 
Planck equation. Computer Methods in Applied Mechanics and Engineering, 194:1513-1526, 2005. 

[48] I. Babuska, J. T. Oden, and J. K. Lee. Mixed-hybrid finite element approximations of second-order elliptic 
boundary-value problems. Computer Methods in Applied Mechanics and Engineering, 11:175-206, 1977. 

[49] J. T. Oden and O.-P. Jacquotte. Stability of some mixed finite element methods for Stokesian flows. Computer 
Methods in Applied Mechanics and Engineering, 43:231-247, 1984. 

[50] T. J. R. Hughes, L. Franca, and M. Balestra. A new finite element formulation for computational fluid dynamics: 
V. Circumventing the Babuska^-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem 
accommodating equal-order interpolations. Computer Methods in Applied Mechanics and Engineering, 59:85-99, 
1986. 

[51] J. Douglilas and J. Wang. An absolutely stabilized flnite element method for the Stokes problem. Mathematics 

of Computation, 52:495 508, 1989. 
[52] L. P. Franca, S. L. Frey, and T. J. R. Hughes. Stabilized finite element methods: I. Application to the advective- 

diffusive model. Computer Methods in Applied Mechanics and Engineering, 95:253-276, 1992. 
[53] L. P. Franca and A. Russo. Unlocking with residual-free bubbles. Computer Methods in Applied Mechanics and 

Engineering, 142:361-364, 1997. 
[54] E. Burman and P. Hansbo. Stabilized Crouzeix-Raviart Elements for the Darcy-Stokes Problem. Numerical 

Methods for Partial Differential Equations, 21:986-997, 2005. 
[55] M. R. Correa and A. F. D. Loula. Unconditionally stable mixed finite element methods for Darcy fiow. Computer 

Methods in Applied Mechanics and Engineering, 197:1525-1540, 2008. 
[56] A. Quarternio and A. Valli. Numerical Approximation of Partial Differential Equations. Springer series in com- 
putational mathematics. Springer- Verlag, New York, USA, 1997. 
[57] H.-G. Roos, M. Stynes, and L. Tobiska. Numerical Methods for Singularly Perturbed Differential Equations. 

Springer- Verlag, New York, USA, 1996. 
[58] Z. Chen, G. Huan, and Y. Ma. Computational Methods for Multiphase Flows in Porous Media. Society for 

Industrial and Applied Mathematics, Philadelphia, USA, 2006. 
[59] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. John Wiley & Sons, Inc., Chichester, UK, 

2003. 

[60] R. E. Ewing and R. F. Heinemann. Convergence analysis of an approximation of miscible displacement in 

porous media by mixed finite elements and a modified method of characteristics. Computer Methods in Applied 

Mechanics and Engineering, 47:73-92, 1984. 
[61] G. Chavent, G. Cohen, and J. Jaffre. Discontinuous upwinding and mixed finite element for two-phase fiows in 

reservoir simulation. Computer Methods in Applied Mechanics and Engineering, 47:93-118, 1984. 
[62] R. E. Ewing and R. F. Heinemann. Mixed finite element approximation of phase velocities in compositional 

reservoir simulation. Computer Methods in Applied Mechanics and Engineering, 47:161-175, 1984. 



18 



[63] L. J. Durlofsky. Accuracy of mixed and control volume finite element approximations to Darcy velocity and 

related quantities. Water Resources Research, 30:965-973, 1994. 

[64] T. Arbogast, M. F. Whcclcr, and I. Yotov. Mixed finite elements for elliptic problems with tensor coefficients as 
cell-centered finite differences. SIAM Journal on Numerical Analysis, 34:828-852, 1997. 

[65] L. Bergamaschi, S. Mantica, and G. Manzini. A mixed finite element-volume formulation of the black oil model. 
SIAM Journal of Scientific Computing, 20:970-997, 1998. 

[66] F. Brezzi, T. J. R. Hughes, L. D. Marini, and A. Masud. Mixed discontinuous Galerkin method for Darcy flow. 
SIAM Journal of Scientific Computing, 22:119-145, 2005. 

[67] E. Burman and P. Hansbo. A unified stabilized method for Stokes' and Darcy's equations. Journal of Computa- 
tional and Applied Mathematics, 198:35-51, 2007. 

[68] J. Bear. Hydraulics of Groundwater. McGraw-Hill, New York, USA, 1979. 

[69] P. Notz, S. R. Subia, M. Hopkins, H. K. Moffat, and D. R. Noble. Aria 1.5 User Manual. Technical report, Sandia 

Report SAND2007-2734, Sandia National Laboratories, 2007. 
[70] K. Devine, E. Boman, R. Heaphy, B. Hendrickson, and C. Vaughan. Zoltan data management services for parallel 

dynamic applications. Computing in Science and Engineering, 4:90-97, 2002. 
[71] A. Ebigbo, H. Class, and R. Helmig. CO2 leakage through an abandoned well: problem-oriented benchmarks. 

Computers & Geoscience, 11:103-115, 2007. 
[72] J. M. Nordbotten, M. A. Celia, S. Bachu, and H. K. Dahle. Semianalytical solution for CO2 leakage through an 

abandoned well. Environmental Science and Technology, 39:602-611, 2005. 



19 



pi 



p2 



1.0 



Figure 1 . A pictorial description of the one-dimensional problem. 
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Figure 2. One-dimensional problem: In the top figure Barus formula is employed, 
and in the bottom figure linear variation of viscosity with respect to pressure is 
employed. Pressure is plotted along x for various values of /?. In this numerical 
example we have taken ao = 1, P{x = 0) = 200, p{x = 1) = 1, and exoL = 10^^". 
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Figure 3. One-dimensional problem: Velocity is plotted along x for various values 
of /3. In this numerical example we have taken ckq = 1, /3 = 0.02, p{x = 0) = 200, 
p{x = !) = !, and exoL = 10"^°. 
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Figure 4. One-dimensional problem: This figure shows the number of iterations 
taken by the numerical algorithm for various values of /3, which is an indicator of 
the strength of the nonlinearity. 
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Figure 5 . One-dimensional problem: This figure shows the variation of the residual 
with respect to iteration number for various values of /3. Note that the y-axis is 
natural logarithm of the residual. 
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Figure 6. Two-dimensional constant flow problem: Four-node quadrilateral (left) 
and three-node triangular (right) meshes used in numerical simulations. 
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Figure 7. Two-dimensional constant flow problem: Barus formula is employed 
with oq = 1, and the tolerance is taken to be etol = 10^^^. For the left figure 
we have used f3 = 0.1 and for the right figure we have used (3 = 0.4. Four-node 
quadrilateral mesh is used for top figures, and three-node triangular mesh is used 
for bottom figures. The numerical results matched well with the analytical solution. 
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Figure 8. A pictorial description of the quarter five-spot problem. 
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Figure 9. Five-spot problem: Pressure contours using four-node quadrilateral 
(top) and three-node triangular (bottom) elements. We have used 21 nodes along 
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each side of the computational domain. 
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Figure 10. Five-spot problem: Variation of Wp'y^ ~ Ph ^^W (which is based on the 
2-norm of the nodal values of the pressure) with respect to iteration number using 
four-node quadrilateral (denoted by Q4) and three-node triangular (denoted by T3) 
elements. In this numerical simulation we have used etol = 10"^*^. 
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Figure 11. Five-spot problem: This figure shows the number of iterations taken 
by the numerical algorithm for various values of (3, and for both Q4 and T3 meshes. 
Note that /3 is an indicator of the strength of the nonlinearity. 
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Figure 12. Five-spot problem: This figure shows the variation of the residual with 
respect to iteration number for various values of /3, and for both Q4 and T3 meshes. 
Note that the y-axis is natural logarithm of the residual. 
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Figure 13. Checkerboard problem: A pictorial description. In regions I and IV, 
we have taken ao = 1-0; and in regions II and III, we have taken ao = 0.001. 
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and B = are employed. 
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Figure 15. Checkerboard problem: Variation of — 



(which is based on 

the 2-norm of the nodal values of the pressure) with respect to iteration number 
using four-node quadrilateral (denoted by Q4) and three-node triangular (denoted 
by T3) elements. In this numerical simulation we have used exoL = 10^^, and 
P = 0.3. 
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Figure 16. Numerical /i-convergence studies: The figure shows the rate of conver- 
gence in L2-norm for four-node quadrilateral (Q4) and three-node triangular (T3) 
elements. 
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Figure 17. Numerical /i- convergence studies: Pressure contours using 21 x 21 four- 
node quadrilateral mesh. We have employed equal-order interpolation for the ve- 
locity and pressure, and there are no spurious oscillations in the pressure. The 
analytical solution for pressure is given by p{x,y) = x'^y'^, and the numerical solu- 
tion matched well with the analytical solution. 
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Figure 18. Three-dimensional constant flow: This figure shows the contours of 
pressure, and the mesh is also shown in the top figure. In the bottom figure we have 
shown the pressure on the x = 1 and x = 4 planes. In this numerical example we 
have employed Barus formula with oq = 1 and (3 = 0.1. 
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Figure 19. Regions with different permeability: This figure shows the computa- 
tional domain and the location of production well. 




Figure 20. Regions with different permeability: Three-dimensional finite element 
mesh using eight-node (linear) brick elements. 
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Figure 21. Regions with different permeability: This figure sliows the velocity 
streamlines for /3 = 0. (Only the streamlines for /3 = are plotted because the 
streamlines for all other values are similar.) 
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Figure 22. Regions with different permeability: This figure shows the production 
rate at the well opening for various values of j3. 
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Figure 23. Regions with different permeability: This figure shows the ratio of the 
total production emanating from regions A and B with respect to (3. 
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Figure 24. CO2 leakage through an abandoned well: A pictorial description of the 
problem showing the cross-sectional view along the plane containing the injection 
and abandoned wells. 
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Figure 25. CO2 leakage through an abandoned wclh This figure shows the com- 
putational mesh used in the numerical simulation. The top-left subfigurc shows the 
mesh at the bottom of the domain, and the top-right subfigure shows the mesh at 
the top of the domain. The bottom subfigure shows the slice through the center 
near the injection and abandoned wells. There are over 1.14 million unknowns in 
this test problem. 
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Figure 26. CO2 leakage through an abandoned well: This figure shows contours of 
the pressure for /3 = (top) and /3 = 1 (bottom). It is evident from the figure that 
the modified Darcy model predicts higher pressures and higher pressure gradients 
than of Darcy model. Half of the domain has been removed to show the detail. 
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Figure 27. CO2 leakage through an abandoned weh: This figure shows contours 
of the magnitude of the velocity for ^ = (top) and /3 = 1 (bottom). It is evident 
from the figure that Darcy model predicts higher velocity in the abandoned well 
than the prediction made by the modified Darcy model. Half of the domain has 
been removed to show the detail. 
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Figure 28. CO2 leakage through an abandoned well: This figure shows the ratio 
of leakage rate to injection rate versus j5. Darcy model over-predicts the leakage 
rate than the modified Darcy model, and the ratio between the leakage rate and 
the injection rate decreases monotonically with increase in /3. This discovery will be 
crucial in designing the seal in a geological carbon-dioxide sequestration geosystem. 
A design of the seal based on Darcy model will be too conservation, which may 
increase the expense of the geological carbon-dioxide sequestration project. 
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